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Abstract. Evolutionary synthesis models are a fundamental tool to interpret the properties of observed stellar 
systems. In order to achieve a meaningful comparison between models and real data, it is necessary to calibrate the 
models themselves, i.e. to evaluate the dispersion due to the discreteness of star formation as well as the possible 
C^h' model errors. In this paper we show that linear interpolations in the log M — log tk plane, that are customary in the 

evaluation of isochrones in evolutionary synthesis codes, produce unphysical results. We also show that some of the 
5_l ■ methods used in the calculation of time-integrated quantities (kinetic energy, and total ejected masses of different 

elements) may produce unrealistic results. We propose alternative solutions to solve both problems. Moreover, we 
have quantified the expected dispersion of these quantities due to stochastic effects in stellar populations. As a 
particular result, we show that the dispersion in the 14 N/ 12 C ratio increases with time. 
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1. Introduction Nevertheless, several aspects of synthesis models are 

still largely perfectible. An analysis of the way in which 



Since their early introduction ( [Tinsley 1980| ), evolution- homology relations of massive stars (implicitly used in 



ary synthesis models have evolved to increase our un- evolutionary tracks interpolations and isochrone computa- 
derstanding of the evolution of stellar populations in the tions) should be modified according to the assumed mass- 
Universe. The improvement of the observational capabili- loss rate, and the inclusion of such new relations into the 
ties has forced the model developers to include more realis- models, still need to be performed. Additionally, the math- 
tic physical ingredients in the models (atmosphere models, ematical approximations used to estimate the lifetimes of 
grids of tracks covering all the evolutionary phases, etc..) massive stars must be carried out carefully, otherwise they 
to interpret the new data. Thus, evolutionary synthesis could produce unphysical results. Finally, the dispersion 
models have become a useful tool to understand the prop- m the model output parameters due to the discreteness of 
erties of observed stellar systems and to test the validity stellar populations has been evaluated only in a few cases, 
of different evolutionary tracks. The improvements of syn- 
thesis models for star-forming regions comprise mainly the Th is work is the third paper of an on-going series 
inclusion of new physical inputs and the extension of the whose the objective to study the oversimplifications and 
output results to more observables. thc possible biases of the evolutionary synthesis models 

for starburst regions, and to assess the confidence limits 

of their outputs. The g lobal structure of the project is 

Send offprint requests to: Miguel Cervirio the following: Paper I flCcrvino ct al. 2000b] ) has been 
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tions. In Paper II (Cervino et al. 2001b) we have investi- 
gated the Poissonian dispersion due to finite populations 
in non-time-integrated observables, its quantitative eval- 
uation and implementation on codes with an analytical 
approximation of the Initial Mass Function (IMF), as well 
as its relation with the dispersion in the output results of 
Monte Carlo simulations. In this paper we present an an- 
alytical approximation to the Supernova rate (SNr) calcu- 
lations in starburst galaxies and the problems related with 
its determination in evolutionary synthesis codes. We also 
study the influence of the interpolations in time- integrated 
quantities (kinetic energy, Eki n , and total ejected masses 
of different elements, y z ) and propose a more precise inter- 
polation technique in order to avoid the unphysical results 
obtained by the previous models. Some improvements of 
the interpolation techniques of evolutionary tracks used in 
synthesis models will be presented in Paper IV (Cervino 
2001 in preparation). We will complete the series with a 
global study of the expected dispersion as a function of dif- 
ferent star-formation laws (continuous and extended star 
formation) . 

In section we present the evolutionary synthesis 
code used. In Section || we present an analytical esti- 
mate of the SNr and how it is computed in evolution- 
ary synthesis models. In Section we show how the re- 
leased kinetic energy and the ejected masses are com- 
puted, and we estimate their Poissonian dispersions and 
bias due to different computation techniques. Finally we 
draw our conclusions in Section pi All the results of this 
paper are available in tabular form in our web server at 



http : //www . laef f . esa . es/users/mcs/ 



2. The evolutionary synthesis model 

Since evolutionary synthesis calculations rely on the prop- 
erties of stars with far more mass values than available 
from stellar evolution calculations, one possible source of 
error is the method used to interpolate between the avail- 
able stellar tracks. 

A few works deal with this problem analytically, either 
proposing analytical formulations for some phases of the 
stellar evolution ( Tout et al. 199q), or using an ana lytical 
population synthesis code (Pluschkc ct al. 2001). One 
advantage of analytical formulations is that the functional 
dependence of the output quantities may be obtained. 

In order to understand and quantify the errors intro- 
duced by synthesis codes for star forming regions, we sum- 
marize the main characteristics of synthesis models with 
non-analytical formulations and the methods used to com- 
pute several parameters. 

1. Most synthesis models interpolate in tables of evolu- 
tionary tracks. Such tables are discrete in their mass 
and time entriesn. When one wants to obtain the lu- 
minosity L, the effective temperature T c ff, or other 



properties, including quantities which change abruptly 
along the evolutionary sequence (e.g. surface abun- 
dances), homology relations are assumed to describe 
with sufficient accuracy the dependence of these quan- 
tities with the initial stellar mass. Therefore, interpo- 
lations in the log M — log Ak plane, where M is the ini- 
tial mass and Ak is a generic stellar property at a given 
evolutionary stage, are usually performed. Additional 
interpolations in the log M — logtfc plane are performed 
to obtain isochrones, and their validity will be exam- 
ined in this paper. Whatever the approach is (fully an- 
alytical or table interpolation), continuity of the stellar 
properties is assumed. The problem of discontinuities 
in evolutionary tracks will be discussed in Paper IV. 
To calculate the integrated properties of the stellar 
population (e.g., the luminosity in a given band) a nu- 
merical integration over the IMF-weighted isochrones 
is always needed. Two main approaches are used: ei- 
ther the IMF is binned into a grid of N initial masses, 
Mi, and then, all stars belonging to the same mass bin 
are assumed to have exactly the same properties, or 
the IMF is sampled with Monte Carlo simulations, and 
then each individual star is evolved and the isochrone 
integration is performed adding the evolved stars. 
To avoid a bias due to the choice of the mass bin, 
a dynamical mass grid can be used within the first 
Schaerer and Vacca 199q, or 



method (Meynet 1995 



199£) 



Starburst99 Leitherer et al 

age, the differences in L and T c ff 



: at each computed 
between two stars 



1 In fact they are discrete in the evolutionary sequence, i.e. 
each point in the table are representative of a given evolution- 
ary stage. 



of initial masses Mi and Mj+i, respectively, are con- 
strained to be lower than a given resolution AL and 
AT c ff. Note that the resulting mass grid will be dif- 
ferent at different ages: the total number of bins N 
and the Mi values vary from one computed age to an- 
other (i.e., it is a dynamical mass grid). Such method 
assures that the H-R diagram is mapped in a contin- 
uous way and all the relevant evolutionary phases for 
the given age (i.e. the isochrones) are included in the 
computations. 

In the case of Monte Carlo simulations, either a high 
number of simulations or a high number of stars in 
each individual simulation are needed to produce a well 
mapped isochrone. As a first order estimation, one sim- 
ulation with 10 5 stars in the mass range 2 - 120 M Q 
is required (as used in Mas-Hesse and Kunth 1991) 
to obtain Ultraviolet and optical luminosities similar 
to those of analytical-IMF models. However, the dis- 
persion of Monte Carlo simulations depends on the 
considered observable. Monte Carlo simulations have 
the advantage of allowing the straightforward compu- 
tation of the standard deviations and the confidence 
levels due to the discreteness of the stellar popula- 
tion, provided the number of simulations performed 
is high enough. An additional advantage is that Monte 
Carlo simulations take into account fast evolutionary 
phases that may be lost in analytical simulations. The 
required number of simulations needed to obtain a sat- 
isfactory estimate of the dispersion can be estimated 
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comparing the mean values of the outputs with the 
results of analytical codes as far as they must to co- 
incide (at least for observables not related with fast 
evolutionary phases). 

In both cases, caution must be taken if the evolution- 
ary tracks present discontinuities (this subject is ad- 
dressed in Paper iv). 
3. The final element of synthesis codes is the age reso- 
lution (or time step) used to compute the integrated 
properties of simple stellar populations, i.e. instanta- 
neous bursts. A sufficient temporal resolution, depend- 
ing a priori on the observable, is needed to assure ac- 
curate convolutions over other arbitrary star-formation 
histories (e.g., the case of constant star formation, sub- 
ject addressed in Paper v). 

For this study we have used the updated version of 



the evolutionary synthesis code presented in Mas-Hesse 
and Kunth (1991[ ); pervino and Mas-Hesse (1994J ) 



The 



updated code includes: 



1. The full set of non-rotating Geneva evolutionary tracks 
including standard (Schaerer et al. 1992, and refer- 



ences therein) and enhanced mass- loss rates (Meynet 



et al. 1994) 



2. The metallicity-dependent atmosphere models for nor 
mal stars from Kurucz (1991), the line blanketed non 
LTE model atmospheres for O stars (CoStar, |Schaerer 



fc De Koter 1997|) and the atmosphere models for 



Wolf-Rayet (WR) stars from |Schmutz et al. (1992| ). 
A numerical isochrone integration using a modified dy- 
namical mass bin, now included in the Dec. 2000 re- 
lease of Starburst99 (Leitherer et al. 199E). We have 



also kept the original Monte Carlo formulation. 

Parabolic interpolations in the log M — log tk plane for 

the isochrones computation^ (see below). 

The computation of the dispersion of all quantities for 

comparisons with low-mass stellar systems (see Paper 

ii). 



study of the IMF has been broadly covered in the astro- 
nomical literature (see the volume of Gilmore ct al. 1998) 
We define the IMF as: 



$(M) 



dN 
dM 



AM' 



(1) 



where a is the IMF slope, A is a normalization factor and 
M the initial mass. This function gives us the probability 
of forming a number of stars in a given initial mass range. 
The widely used Salpctcr's IMF slope corresponds to a = 
2.35 with this definition. The number of stars N sta r(t) 
present in a system at a given time from the onset of star 
formation within an initial mass range, is obtained by the 
convolution of $(M) with the Star Formation Rate law, 



ft fM(t) 
N s tar = / / $(M)#(t - t')dMdt', 

JO JM bw 



(2) 



where M(t) is the initial mass of the star that ends its 
evolution at time t. 

Throughout this work we assume an Instantaneous 
Burst (IB) of star formation, ^(t) = S(t) where S(t) is 
Dirac's delta function, so that: 



f M(t) 
N star = / <S>(M)dM. 

iM,„ w 



(3) 



The number of stars that will end their evolution in 
the systemR N$n, in a time interval [ii,i2] is 

N SN = / ${M)dM. (4) 

JM(ti) 

For mathematical convenience M(t) can be approxi- 
mated by 



M{t) = Br T , 



(5) 



with 7 > 0, thus assuming implicitly a linear relation 
between logM and \ogt. Table l] shows the values of 7 
and log B for several mass ranges from the Schallcr ct 



The calculations in the present paper are done with the 



solar metallicity tracks of 3challcr et al. (1992), adopting ten as: 



al. (1992 ) solar metallicity tracks using a linear log M — 
log t approximation (but see below) . Eq. |l| can be rewrit- 



standard mass-loss rates, a Salpeter IMF over the range 
from 2 to 120 M , and an instantaneous burst. 



3. The supernova rate 

The supernova rate, as any other output of an evolution- 
ary synthesis code, depends on the assumed IMF. The 



N SN 



$[M(t)} 



dM 



dt 



dt 



(0) 



Using Eq. Wand EL we obtain the SNr, i.e. the number 
of SN in a time interval: 

In the general case of a function M(t), we obtain the 
exact value of the SNr, i.e. the number of SN in a time 
interval: 



2 During the work on the present models an error, originating 
from the change from track to isochrone synthesis, was found 
in the calculation of t he SNr and some der ived quantities in 
the Starburst99 code (Leitherer et al. 1999). This resulted in 



SNr(t) 



dN SN 
dt 



AM(ty 



dM 



dt 



(7) 



a strongly non-monotonic SNr, partly increasing with time, 
in contrast with the expectations. The method described here 
has now been included in the December 2000 release of the 
Starburst99 code. 



3 This treatment is general for all stars described by the 
function M(t). Since we restrict ourselves to times shorter than 
20 Myr, all the stars considered will end their lives either with 
a SN explosion or with the formation of a Black Hole. Here we 
do not distinguish between these cases. 
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Mass range (Mq) 


7 (Eq. PI) 


log 


B(Eq. 


D 





120 


85 


4.51 




31.31 




5.09 


85 


60 


1.88 




14.15 




1.54 


60 


40 


1.95 




14.61 




1.63 


40 


25 


1.21 




9.68 




0.63 


25 


20 


0.94 




7.81 




0.26 


20 


15 


0.81 




6.96 




0.10 


15 


12 


0.68 




6.04 




-0.08 


12 


9 


0.57 




5.22 




-0.23 


9 


7 


0.50 




4.68 




-0.33 



Table 1. Values of 7 and log-B defined in Eq. [5] and the 
SNr slope (assuming a Salpeter IMF slope) for different 
mass ranges from the Schaller et al. (1992 ) solar metal- 
licity tracks. 
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Using Eq. (i.e. a linear interpolation in log M — log t) 
one obtains: 



SNr(t)=AB- a+1 jt 13 , 



(8) 



where (3 = 7 a — 7 — 1. For the Salpeter IMF Eq. H shows 
that the SNr is a decreasing function of age. This expres- 
sion is also useful to verify the proper calculation of the 
SNr in evolutionary synthesis models that use linear in- 
terpolations in logM — logi. 

3.1. Implementation in evolutionary synthesis codes 

To calculate the SNr in evolutionary synthesis codes, we 
compute the population at some given age, tj. The basic 
idea is to know how many stars have ended their evolution 
between the previous computed age, tj—i, and the current 
one, tj. Then, the SNr obtained is the mean value of the 
SNr for the used time interval. At the age tj this is given 
by: 



SNr(tj 






NsNJt,) 



ti - 1 



3-1 



U -t 



j'-i 



(9) 



where Wi is the normalized number of stars of initial masqj 
Mi, and a,- is defined as 



s if i > i{tj) or i < i(tj~i) 
~~ *• 1 if i(tj-i) < i < i{tj) 



(10) 



where i(tj) is the index in the binned IMF grid corre- 
sponding to a mass M(tj). The indexes are given by the 
function M(t) described above. 

The resulting SNr using a linear or parabolic inter- 
polation of M(t) is shown in Fig. pi Whereas the SNr 
using the linear interpolation (cf. Eq. ph exhibits discon- 
tinuities corresponding to the discreteness of the stellar 
tracks (cf. Table m , the parabolic interpolation presents a 
much smoother behavior. When linear interpolations are 
used, the resulting lifetimes of stars at both sides of a 
given tabulated stellar track are lager than the lifetimes 

4 Note that for codes that use a dynamical mass binning, Wi 
is in fact Wi(t) 



og Age [yr] 

Fig. 1. SNr using different interpolation techniques. The 
solid line corresponds to a linear interpolation in the 
logAf — logi plane. The short-dashed line corresponds to 
a parabolic interpolation. 




Fig. 2. logA/- 
[Bchallcr et al. 



5 10 20 50 100 

Mass [M ] 

log t e plane for two different set of tracks: 
(1992|) at solar metalli city and standard 



mass-loss rate and Meynet et al. (1994) at twice the solar 
metallicity and twice the mass-loss rate 



obtained with parabolic interpolations. This produces a 
lower SNr when the stars with masses corresponding to 
the tabulated track had just exploded and an accumula- 
tion of SN events just before the lifetime corresponding to 
the following tabulated star. 

The situation is also illustrated in Fig g where the 
logM — logi plane is shown for the tracks used and 
the solar-mctallicity and twice solar-metallicity tracks of 
Meynet et al. (1994) (for masses lower than 12 Mq we 
have completed the table with schaerer et al. 1992). 



A general deviation from linearity is present for massive 
stars, such deviation is more extreme in the case of twice 
solar metallicity tracks with twice mass-loss rates. 

It is interesting to remark that whatever the interpo- 
lation technique is, some wiggles are seen at the beginning 
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of the evolution of the SNr. Whereas the abrupt disconti- 
nuities at the ages that correspond to the lifetime of the 
60 and 40 M Q stars are due to the interpolation tecniquc, 
the wiggles themselves are due to the particular behavior 
of lifetime of the WR stars in the set of tracks used. A de- 
tailed analysis of the lifetimes of the tabulated stars that 
reach the WR phase reveals a nonmonotonic behavior of 
the slope of the lifetime itself. This behaviour, convolved 
with the IMF slope, produces the wiggles present in the 
figure; i.e., the presence of the wiggles arise naturally given 
the set of tracks used, whereas the exact ages where the 
wiggles appear are dependent on the interpolation tech- 
nique. 

The output of a synthesis code should not depend on 
the specific masses tabulated in the evolutionary tracks 
(unless the tracks follow a discontinuity of the stellar evo- 
lution). But the behavior of the l^jrl relation shows such 
dependence if a linear interpolation in the log M — logt 
plane is used. Moreover, the linear behavior assumed in 
the log M — log t plane is not real at all for some set of 
tracks. The misbehavior of linear interpolations may be 
due to the effects of mass loss and overshooting in mas- 
sive stars and it may also be present in the new generation 
of tracks with rotation. 

Finally, we remark again that we have focused on the 
ages of the SNe explosions, but a similar situation ex- 
ists at other evolutionary phases. We want to stress that, 
even if the parabolic interpolation used here seem to pro- 
duce more realistic results, a correct interpolation tech- 
nique (based on physical principles) does not yet exist, 
and a more careful study is necessary on this subject. 
The parabolic interpolation subroutines developed here 



are available at tittp : //www. laeff .esa.es/users/mcs/ 



3.2. Estimate of the SNr dispersion obtained from an 
evolutionary code 

The SNr is not a direct observable. However it enters in 
the calculation of other observables like the non-thermal 



radio flux (Mas-Hesse and Kunth 1991) or the ejection 



of elements into the ISM. So, the knowledge of the SNr 
dispersion due to the discreteness of the stellar population 
is needed to obtain the expected dispersion in the observed 
properties of real systems. In the following paragraph we 
summarize how to calculate such a dispersion. We refer to 
Buzzoni (1989 ) and Paper n for further details. 



The IMF gives the probability, Wi, of finding a num- 
ber of stars within a given mass range at t = 0. If we 
assume that each Wi follows a Poissonian distribution, the 
variance of each Wi, of , is equal to the mean value of the 
distribution, Wi. Let us assume now that each star has 
a property a*, so that the contribution to the integrated 
property A of the star of the same mass is given by WiQi 



observable A is the sum of all the variances. The relative 
dispersion is: 



A 






2)1/2 



1 



S*=i w * a * \/N c ft{A) 



(11) 



where the last term gives u s the definition of N g (A) de- 
scribed by Buzzoni (1989| ). Note that N c f[(A) is normal- 
ized to the total mass, so N B g(A) gives directly the relative 
dispersion for any total mass transformed into stars. 

Paper II shows that the dispersion obtained from Eq. 
O is equal to the dispersion of Monte Carlo simulations. 
Let us stress that such dispersion is present in Nature 
(star populations are always discrete and finite) and it is 
not an evaluation of the errors of the synthesis models, i.e. 
the dispersion is also an observable. This intrinsic disper- 
sion must be taken into account, before establishing any 
conclusion, when fitting observed quantities to model out- 
puts. Finally, the evaluation of the dispersion depends on 
the interpolation techniques used, so a correct interpola- 
tion technique is also required to fit this observed property 
of Nature. 

In the case of the NsN{tj), a,; is defined by Eq. |l(J, and 
the relative dispersion in NsN{tj) is: 



1 



a NsN(tj) 

N SN (t 3 ) ~ y/Nsffitj) 



(12) 



with a variance aja? 



Wia: 



? The total variance of the 



In this particular case (a, = or 1) the mean value 
and the variance coincide as it is the case in Poissonian 
distributions. The obtained N$n is a mean value over the 
time step used, and the obtained dispersion shows the 
variation about that mean value. The dispersion, however, 
depends on how the mean value is computed, i.e. depends 
on the time step used. 

The dispersion on the mean SNr is obtained dividing 
the variance <J% S (tj) by the time-step, i.e. N c g(SNr) — 
SNr(tj). Figure || shows the 90% confidence limit for dif- 
ferent Monte Carlo simulations. We have used 500 simula- 
tions of clusters with a total mass transformed into stars of 
10 4 M Q , 200 simulations of clusters with 10 5 M Q , and 100 
simulations of clusters with 10 6 M Q . A time step At = 0.1 
Myr has been used in these simulations. 

The dispersion obtained from the Monte Carlo sim- 
ulations is given by OMon = SNr(t) x At. This time- 
dependence in the dispersion must be taken into account 
for quantities related to the SNr. Note that the obtained 
dispersion is correct if the SNr is defined in units of num- 
ber of SN each 10 5 years instead the usual units of SN per 
year. 

As an example, it is assumed that our Galaxy have a 
mean SNr about 3 SNe per century which means, assum- 
ing a Poissonian distribution for the SNr, asNr ~ 2 SN 
per century and a relative dispersion of 0.6. If the SNr 
is defined as 3 10 4 SN per Myr, the corresponding asNr 
becomes 173 SN per Myr and the relative dispersion is 
0.006. 
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10° Mo, N=100 



10° Mo, N = 200 



10* Mo, N = 500 



Fig. 3. 90% Confidence limit from Monte Carlo simula- 
tions for the SNr of clusters with different amount of total 
mass transformed into stars since the beginning of the 
burst of star formation with a time step of 0.1 Myr 



4. Kinetic energy and ejected masses 

The kinetic energy and ejected masses have two compo- 
nents: stellar winds and SNe. As we have pointed out be- 
fore, all the outputs are affected by the way the interpola- 
tions are performed in the logM — log £& plane. In the fol- 
lowing we use parabolic interpolations in this plane. Note 
that such interpolations produce a bit lower life-time that 
linear ones, which means a lower amount of kinetic en- 
ergy and integrated ejected masses. We now study these 
properties in more detail. 

We assume a typical value of 10 51 erg SNe -1 for the ki- 
netic energy released by a SN. The kinetic power, Pkin{t), 
is the product of the typical energy released by a SN 
multiplied by the SNr(t). In the case of ejected masses, 
we need to use a relation between the SN ejected masses 
and the mass of the exploding star. Illustrative examples 



al. (2000a) 



can be found in Portinari et al. (1998) and Cervino ct 



4.1. Stellar wind components 

The kinetic energy and ejected masses also include a con- 
tribution from stellar winds. The stellar mass-loss rate is 
the key parameter needed to compute the kinetic energy 
and the ejected masses from stars before the end of their 



evolution. The total kinetic power is the sum over the con- 
tribution of individual stars, Pi{i) at age t: 



N 
Pkin{t) = y^^WjPi(t) 



1 



N 

2^ 



w l fhi{t)v 2 00 



i(t), 



(13) 



where m, (t) is obtained from interpolations of the tracks 
and Voo t i(t) is derived from the interpolated luminosity, 
the effective temperature, the mass rrn, and the metallicity 
of the star of initial mass Mi at the given age, following 



Lcithcrcr and Hcckman (1995) 



Similarly, the instantaneous ejected mass of an element 
z, y' z (t), can be computed as: 



Vz® 



N 



N 



WiV. 



i,i(t) = /]Wirhi(t)Zi(t), 



(14) 



where Z$ is the mass fraction of the surface abundance of 
the element z for a star of initial mass Mj at age t. In 
general, in evolutionary synthesis codes, interpolations in 
L, T e ff, rh and m are performed in the logM — log Ak 
plane, and those in Z in the logM — Zk plane, where the 
subindex k refers to a given evolutionary stage. In the 
case of chemical evolution models the interpolations in Z 
vary for different authors, from linear interpolations in the 
M — Z plane ( Fcrrini ct al. 1994 ) to the use of splines 
( ICarigi 2000| ). 

For the computation of the corresponding time inte- 
grated quantities — the cumulative kinetic energy Ekm 
and the total yield y z — an additional sum over time is 
needed. Two different approximations may be followed: 

— Method (a) The results from previous computed ages 
are used to compute the sum: Each value of Pkin(tj) or 
y z (tj) can be considered constant over the time interval 
between ij_i and tj, where the index j defines the age 
array used for the code output, and the Ekm and y z 
are obtained adding up such contributions. Note that 
this method depends on the age array used and will 
lose evolutionary stages where the characteristic time 
is shorter than the time step used. 

— Method (b): A similar treatment to those of chemical 
evolution models is used. An additional table for each 
evolutionary state, k, with the integrated amount of 
kinetic energy or chemical abundances from t = to 
each tabulated age, tk is computed. The final point 
of the table is the integrated energy and ejected mass 
resulting from the action of stellar winds all along the 
evolution of the star. 

This method has two possible implementations: 
— Method (b.l): The kinetic energy or chemical abun- 
dances tables are obtained using the mass-loss, i.e. 
for the chemical abundances case: 

k' 



y Z ,i( t k<) = ^2 "*<(**) Z ii t k) [U 



ti.k- 



k-l, 



(15) 



fc=i 



where the subindex k refers to the tabulated points 
in the tracks for the star of initial mass M, . 
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Method (b.2): The instantaneous mass is used in- 
stead of the mass loss: 



Vz,i{tw) 



k' 

E 

fc=i 



m,i(tk) - mi(tfc-i)] < Zi(t k ) >, (16) 



where < Zi(tk) > is the mean value of the surface 
abundance at the evolutionary ages t k and ifc-i of 
a star of initial mass M,-. This is the approximation 
we have used here. 

In both cases, the total amount of the ejected element 

z at a given age t becomes: 



Vz(t) 



N 

E 



WiVz,i(t) 



(17) 




This method has the advantage that the output does 
not depend on the age array used in the synthesis code. 

The three methods should converge to the same value, 
but th is can be only achieved if the time step used in 



10 20 

Initial Mass [M Q ] 

Fig. 4. Ratio of the integrated mass-loss during the life- 
time of the star vs. "reconstructed" mass loss by subtrac- 
tion of the mass at the Zero Age Main Sequence (ZAMS) 
and the mass at the end of the evolution for Schallci 



metho ql (a) is the same as the lowest time step used in the 
evolutionary tracks (i.e. a few years, that it is prohibitive 
for realistic computations). Methods (b.l) and (b.2) must 
also converge if 



the mass at the end of the evolution for 
ct al. (1992J ) tracks (solid line) at solar metallicity and 
standard mass-loss rate and [Mcynct ct al. (1994 ) tracks 



at twice solar metallicity and twice mass-loss rate tracks 
(dashed line). 



rhiitk) {t>Lk - U,k-i) = m i(tk) - mi(*fc_i), 



(18) 



which is, however, not found to be true for various sets of 
stellar tracks adopted. 

Let us illustrate the situation defining the ratio R for 
each star as: 



R = 



ZfcLi™fa)fa -*fc-i) 
M - m ke 



(19) 



where the index k e refers to the last tabulated point in the 
track. Such ratio must be equal to 1 if Eq. ^8| is fulfilled. 
The resulting values are shown in Fig. and in Table for 
the Bchaller ct al. (1992) tracks at solar metallicity and 
standard mass-loss rate and Meynet et al. (1994) tracks 



Initial 


J^mAt 


M - m fce 


J^rn-At 


M - m ke 


Mass range 


(b.l) 


(b.2) 


(b.l) 


(b.2) 


(M ) 




Z Q 


2xZ Q 
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12 


0.5 


0.4 


2.3 


2.3 


15 


1.5 


1.4 


3.5 


3.5 


20 


3.6 


3.5 


10.1 


9.8 


25 


9.6 


9.4 


18.6 


19.3 


40 


31.3 


31.9 


33.0 


35.5 


60 


47.6 


52.2 


63.3 


56.2 


85 


68.3 


76.0 


83.8 


82.5 


120 


129 


112 


120 


118 



at twice solar metallicity and twice mass- loss rate. 

Typically, as shown m Fig. and Table Q the ratio 
R between the "reconstituted" integrated mass loss using 
the above methods and the total mass loss given by the 
difference between the tabulated initial and final mass of 
the tracks, is found to vary by up to ~ 10 % for stars 
with non negligible mass loss. Note that for a 120 M Q 
star the integrated mass loss is equal or higher than the 
initial mass! 



4.2. Global evolution and dispersion 

We now study the resulting integrated kinetic energy and 
chemical yields including both stellar winds and SNe. As 
an example, we have focused on the ejected mass of 12 C 
and 14 N/ 12 C ratio. These elements, principally 14 N, in 
a "standard" stellar population are mostly produced by 
the intermediate stars in the Asymptotic Giant Branch 
(AGB) phase and when they form Planetary Nebulae 



Table 2. Values of the mass loss between the ZAMS 
and the end of the evolution, and the integrated mass 
loss during the lifetim e for different mass ranges f rom the 
[Schallcr ct al. (1992 ) solar metallicity tracks and Meynet 



ct al. (1994) tracks at twice solar metallicity and twice 



mass-loss rate. 



(PNe). However AGB stars and PNe appear at later ages 
than those discussed in our work. Therefore, the 14 N and 
12 C produced in the first few Myr come only from mas- 
sive stars through stellar winds and SN explosions. We 
have selected these two elements as illustrative examples, 
to highlight the importance of the WR-wind phase in mas- 
sive stars. 

Figures pL ra and R show these quantities computed 
with method (a) for various time steps At, and using 
method (b.2) with different interpolation techniques. For 
method (a) we use linear interpolations in log M — A for 
the abundances, and log M — log A for E k i n - For method 
(b.2) we compare linear interpolations in the M — A and 
log M— log A planes. The examination of the figures shows 
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Fig. 5. Top: Integrated kinetic energy using three different 
time steps with method (a) and using method (b.2) de- 
scribed in the text with interpolations in the M — E^m 
plane (linear) and logM — log Ekin plane (log - log). 
Bottom: N c g(Eki n ) as defined in Eq. 11 and computed 
by method (b.2) using the M — Ekin plane for interpola- 
tions. The right vertical axis shows the minimum amount 
of gas that needs to be transformed into stars (in the given 
mass range and for the given IMF slope) to give a relative 
dispersion lower than 10% when analytical-IMF models 
are compared with real data. 



the following: i) time steps Ss 0.1 Myr using method 
(a) appear adequate to properly calculate the considered 
quantities; ii) at young ages dominated by stellar winds 
(t < 3 Myr), the use of pretabulated integrated quantities 
(method b) leads to a somewhat larger and more correct 
values for kinetic energy and yields produced by massive 
stars. The origin of this difference is due to rapid varia- 
tions of these quantities along the isochrone, which are not 
well enough resolved by method (a). The differences be- 
tween the different numerical techniques are relevant for 
the younger ages, when the integration of stellar winds 
are the only contribution of the time-integrated quanti- 
ties. Using method (a) our code needs about 30' of CPU 
time in a SunOS spare machine using a time step of 0.1 
Myr form 0.1 to 20 Myr. The CPU time is in this case in- 
versely proportional to the time step, i.e. the code would 
require about 300' of CPU time to obtain very similar re- 
sults once the SN activity is the dominant source, with a 
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Fig. 6. Integrated amount of 12 C ejected using three dif- 
ferent time steps with method (a) and using method 
(b.2) described in the text with interpolations in the 
M — 12 C plane (linear) and log M— log C plane (log - log). 
Bottom: N cS ( 12 C) from method (b.2) using the M -C 
plane for interpolations. The right vertical axis shows the 
minimum amount of gas that needs to be transformed 
into stars (in the given mass range and for the given IMF 
slope) to give a relative dispersion lower than 10% when 
analytical-IMF models are compared with real data. 



time step of 0.01 Myr. On the other hand, the time com- 
putations by method (b.2) are not time-step dependent 
and they produce stable results whatever the time step is. 

In the lower panel of Fig. ||, || and we show the result- 
ing iV c ff for the kinetic energy and ejected masses resulting 
from method (b.2) using interpolations in the M — A plane. 
N c ff can be easily computed for the integrated properties 
if method (b.2) is used. Note, though, that the evaluation 
N e g of the integrated properties from method (a) with a 
dynamical mass binning is quite difficult because we need 
to know the contribution of the same individual popula- 
tion along all the computed ages. 

In the case of the kinetic energy, the relative dispersion 
decreases with age once the first SNe appears in the clus- 
ter. The natural explanation is that more and more stars 
contribute to the released kinetic energy and the statistics 
becomes better and better. 

For the 12 C case, there are three dominant contribu- 
tions at different times: massive non-WR stellar winds 
from to 3 Myr, WR stellar winds from 3 to 5 Myr and 
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the SNe thereafter. The depression in N cS ( 12 C) at 3 Myr 
coincides with the onset of the WR phase. As it has been 
shown in Paper II, the WR phase is strongly influenced 
by the discreteness of the stellar population and shows a 
high dispersion in the related properties. The cumulative 
production of 12 C by SNe decreases again the dispersion 
after 5 Myr. 

The 14 N/ 12 C ratio (Fig. [?], left panels) shows several 
interesting effects. It is characterized by an increasing ratio 
until 3 Myr due to the effects of stellar winds of massive 
stars. WR stars start to appear in the cluster at about 
2 Myr. The evolution of WR stars at solar mctallicity 
follows the sequence of WR stars with N in the envelope 
(WN phase, characterized by a strong mass-loss rate) and 
a posterior WC phase, where C appears at the surface in 
an amount larger than N, with a mass-loss rate dependent 
on the mass of the WR (i.e. decreasing with time). The 
prevalence of massive OB stars and the WN phase last 
until 3 Myr and produce more N than C. Later the action 
of WC winds and SNe increases the C production. 

The top right panel of Fig. [7] shows the evolution of 
iV e ff( 14 N/ 12 C). The relative dispersion (inverse of N e g) 
reaches a minimum value during the WR phase as for 
12 C. N B g increases from the first SN explosion to 5 Myr, 
and it decreases again slowly for more evolved ages. The 
computation of iV e ff( 14 N/ 12 C) must take into account the 
effect of the covariance (i.e. the correlation coefficient, 
p( 14 N/ 12 C)). Due to correlation effects, the dispersion in 
the 14 N/ 12 C ratio is lower than the dispersion on 12 C or 
14 N alone. Note also that the dispersion of the 14 N/ 12 C 
ratio increases as the cluster evolves. 

Finally Fig. shows the resulting values of Ekin and 
14 N/ 12 C ratio from a set of Monte Carlo simulations. 
The figure shows the 90% confidence level for simulations 
where different amounts of gas have been transformed into 
stars. The exact values of the 90% confidence level can not 
be achieved by the numerical formulation proposed, as far 
as the corresponding probability density distributions are 
not known, but it is clear that the dispersion in the simu- 
lations have a similar behavior than the one obtained by 
the computation of N c g. In particular, Monte Carlo sim- 
ulations and the computation of N e g ( 14 N/ 12 C) show that 
the dispersion in the 14 N/ 12 C ratio increases with age. 

5. Conclusions 

We have discussed technical issues (accuracy, impact of 
various interpolation and integration methods) regarding 
the calculation of supernova rates, as well as chemical 
and mechanical yields in evolutionary synthesis models. 
Furthermore we have quantified the expected dispersion 
of these quantities due to stochastic effects in populations 
of various total masses. The main conclusions are the fol- 
lowing: 

1. Linear interpolations in the log M — log tk plane, where 
M is the initial mass and t is the age for a given evo- 
lutive stage, can give unphysical results for the super- 





Fig. 8. 90% Confidence level of Monte Carlo simulations 
for Ekin and 14 N/ 12 C ratio of clusters with different 
amount of total mass transformed into stars since the be- 
ginning of the burst of star formation. The grey scale val- 
ues are the ones defined in Fig. 



nova rate in particular and for the evolutionary phases 
in general 

A parabolic interpolation technique can improve the 
results, but a more general technique, taking into ac- 
count the stellar evolution theory (i.e. the relation of 
the stellar evolution parameters like luminosity and 
effective temperature with time) would be an asset. 
However, as far as stellar theory is not complete, 
parabolic interpolations can be used to produce more 
reasonable results than linear ones. 
The unphysical results produced by the linear interpo- 
lations do not only affect the massive stars (as naively 
expected) , but also the low mass ones (at least down to 
9 M Q ). The effects must be quantified for a wide age 
range, including low mass stars, in the evolutionary 
synthesis models. 

The time-integrated quantities of instantaneous burst 
models (or single stellar populations) depend on the 
time step used for the integrations. For the quantities 
considered here a time step not larger than ~ 0.1 Myr 
is found to be required for the integration of the stel- 
lar winds component, that are relevant at the early 
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Fig. 7. Top left: 14 N/ 12 C ratio obtained from the ejected 14 N and 12 C computed for three different time steps with 
method (a) and using method (b.2) described in the text with interpolations corresponding to Fig. a Bottom left: 
Detail of 14 N/ 12 C ratio obtained at the beginning of the burst of star formation. Top right: 7V e ff( 14 N/ C) (see text). 
The right vertical axis shows the minimum amount of gas that needs to be transformed into stars (in the given mass 
range and for the given IMF slope) to give a relative dispersion lower than 10% when analytical-IMF models are 
compared with real data. Bottom right: Correlation coefficient, p( 14 N/ 12 C). 



phases of the evolution of the cluster. However, the 
best choice in terms of computing time and accuracy 
is the use of predefined tables with the corresponding 
time-integrated quantities. 

In early phases (t Ss 3 Myr) the kinetic energy and the 
ejected elements of star forming regions depend very 
strongly on the stellar winds. In these phases, the most 
accurate outputs are obtained when time-integrated 
tables of evolutionary tracks are used. 
The discreteness of the real stellar populations is ex- 
pected to produce a dispersion in the observed param- 
eters that must be taken into account a priori when 
compared to the outputs of synthesis models. For time- 
integrated quantities, the dispersion is higher at the 
beginning of the star formation episode and becomes 
even more important during the WR phase. The de- 
pendence of the theoretical dispersion on the total stel- 
lar mass has been quantified. 



5. When the correlation between different yields of ele- 
ments is taken into account, the dispersion of the ratio 
of such elements increases with time. The relevance of 
this effect on the observed dispersion of the 14 N/ 12 C 
ratio remains to be evaluated. 
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